Determine temperature rise stressor per spp

knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)

library(raster)
library(sf)
library(tidyverse)
library(here)
library(RColorBrewer)
library(oharac)

source(here('common_fxns.R'))
source('stressor_fxns.R') ### in same dir as this Rmd

1 Summary

Using a map of mean SST over the past five years (2016-2020), identify species-specific stressor based on the species’ thermal envelope (determined from AquaMaps envelopes or IUCN-similar envelopes).

For each unique thermal envelope (i.e., multiple species might share the same absolute max and preferred max):

• Identify the base stressor layer using Mel’s function to convert mean temperature to 0-1 stressor
  • If mean temp exceeds absolute max, stressor = 1
  • If mean temp is below preferred max, stressor = 0
  • Otherwise stressor scales linearly with temperature from preferred max to absolute max
  • Not going to worry about temps below preferred/absolute min temps
• For each species with this thermal envelope, mask the base stressor layer with the species range and save out as .csv

2 Data

• Species range maps from AquaMaps and IUCN
• Species thermal envelopes from AquaMaps and calculated for IUCN
• Temperature rise data from ???

3 Methods

3.1 Read in mean temp map

This is stored on Mazu at /home/shares/ohi/stressors_2021/_dataprep/SST_past_present/final_avg_tmp/sst_avg_2016-2020.tif.

mean_T_rast <- raster('/home/shares/ohi/stressors_2021/_dataprep/SST_past_present/final_avg_tmp/sst_avg_2016-2020.tif')
plot(mean_T_rast, axes = FALSE, main = 'Mean annual T, °C')

# class      : RasterLayer 
# dimensions : 1814, 3617, 6561238  (nrow, ncol, ncell)
# resolution : 10000, 10000  (x, y)
# extent     : -18086282, 18083718, -9069952, 9070048  (xmin, xmax, ymin, ymax)
# crs        : +proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
# source     : sst_avg_2016-2020.tif 
# names      : sst_avg_2016.2020 
# values     : -1.548435, 30.61173  (min, max)

3.2 Determine thermal envelope for species

For AquaMaps read from their envelope data (using get_am_spp_envelopes() helper function); for IUCN, read from results of script 6a_thermal_tolerances_estimate_iucn.Rmd.

Here we are calculating for all AquaMaps and IUCN maps, regardless of whether there is overlap or the number of occur_cells is too low. We can always ignore those maps later.

am_spp_info <- get_am_spp_info()

am_T_env <- get_am_spp_envelopes() %>%
  filter(param == 'temp' & dist %in% c('max', 'pref_max')) %>%
  spread(dist, value) %>%
  left_join(am_spp_info, by = 'am_sid') %>%
  select(am_sid, sciname, T_max_a = max, T_max_p = pref_max)

iucn_spp_info <- read_csv(here('_data/iucn_spp', 'iucn_to_worms_match.csv'))
iucn_maps <- read_csv(here('_data/iucn_spp/spp_marine_maps_2021-3.csv'))
iucn_derived_T_envs <- read_csv(here('_setup/stressors/int/iucn_spp_thermal_env.csv')) %>%
  select(iucn_sid, iucn_tmaxa = T_max_a, iucn_tmaxp = T_max_p)

iucn_T_env <- iucn_spp_info %>%
  filter(iucn_sid %in% iucn_maps$iucn_sid) %>%
  ### first try to join AquaMaps temp profiles
  left_join(am_T_env, by = c('worms_name' = 'sciname')) %>%
  left_join(iucn_derived_T_envs, by = 'iucn_sid') %>%
  mutate(T_max_a = ifelse(is.na(T_max_a), iucn_tmaxa, T_max_a),
         T_max_p = ifelse(is.na(T_max_p), iucn_tmaxp, T_max_p)) %>%
  select(iucn_sid, sciname, T_max_a, T_max_p) %>%
  filter(!is.na(T_max_a) & !is.na(T_max_p))

T_env_all <- bind_rows(am_T_env, iucn_T_env)

hist(T_env_all$T_max_a, main = 'T max, absolute')

hist(T_env_all$T_max_p, main = 'T max, preferred')

hist(T_env_all$T_max_a - T_env_all$T_max_p, main = 'T max, abs - pref')

3.3 Process species pressures

Group species into similar thermal niches. Round T_max_p down to the nearest 0.5°, round T_max_a up to the nearest 0.5°, then group by T_max_a and T_max_p. This floor/ceiling rounding will expand the thermal envelope slightly, resulting in a more conservative estimate of pressure. Grouping like this results in 1385 distinct groups with similar thermal niches.

For each niche group:

• calculate the pressure map using the calc_T_max_prs() function defined below
• convert the resulting pressure map to a data.frame based on cell ID and mean annual temp
• loop over each species in the niche:
  • read in the species distribution as a .csv from the appropriate source (AquaMaps or IUCN)
  • intersect the species distribution with the temperature pressure map for this thermal niche
  • write out the resulting csv of cell_id and pressure, with species info in filename.

3.3.1 Function calc_T_max_prs()

To determine the pressure on a species based on its thermal tolerance and mean annual temperature in a cell:

• If the mean annual temperature meets or exceeds the species’ T_max_a, pressure = 1.
• If the mean annual temperature meets or is lower than the species’ T_max_p, pressure = 0.
• If the mean annual temperature T_max_p < mean annual T < T_max_a, scale linearly from 0 to 1.

### before rounding, 10175 gps
### * after rounding to 0.1° 6191 gps
### * after rounding to 0.25° 3397 gps
### * after rounding to 0.5° 1415 gps
### * rounding pref down and abs up (floor/ceiling), 462 gps
# gp_check <- T_env_all %>%
#   # mutate(T_max_a = round(T_max_a, 1), T_max_p = round(T_max_p, 1)) %>%
#   # mutate(T_max_a = ceiling(4 * T_max_a) / 4, T_max_p = floor(4 * T_max_p) / 4) %>%
#   mutate(T_max_a = ceiling(2 * T_max_a) / 2, T_max_p = floor(2 * T_max_p) / 2) %>%
#   # mutate(T_max_a = ceiling(T_max_a), T_max_p = floor(T_max_p)) %>%
#   group_by(T_max_a, T_max_p) %>%
#   summarize(n = n())

### For synonymous spp with multiple AquaMaps instances, take mean of temps
T_env_rounded <- T_env_all  %>%
  mutate(src = ifelse(is.na(am_sid), 'iucn', 'am'),
         id  = ifelse(src == 'iucn', iucn_sid, str_replace_all(sciname, '[^a-z]+', '_'))) %>%
  group_by(sciname, id, src) %>%
  summarize(T_max_a = mean(T_max_a), T_max_p = mean(T_max_p),
            .groups = 'drop') %>%
  mutate(T_max_a = ceiling(2 * T_max_a) / 2, T_max_p = floor(2 * T_max_p) / 2)


gps_df <- T_env_rounded %>%
  group_by(T_max_a, T_max_p) %>%
  summarize(n_spp = n(), .groups = 'drop') %>%
  arrange(desc(n_spp)) %>%
  mutate(gp = 1:n())

T_env_to_process <- T_env_rounded %>%
  left_join(gps_df, by = c('T_max_a', 'T_max_p')) %>%
  filter(!is.na(sciname)) %>% 
  distinct()

### zxcv <- show_dupes(T_env_to_process, 'sciname')
### qwer <- show_dupes(T_env_to_process, 'id')

calc_T_max_prs <- function(abs, pref, rast = mean_T_rast) {
  #t0 <- Sys.time()
  pressure <- (1/(abs - pref))*(rast - pref)
  pressure[pressure < 0] <- 0
  pressure[pressure > 1] <- 1
  return(pressure)
  #Sys.time()-t0
}

# tmp <- calc_T_max_prs(abs = 30, pref = 5, rast = mean_T_rast)
# plot(tmp)

get_spp_map <- function(map_f) {
  
  if(!file.exists(map_f)) stop('Map does not exist: ', basename(map_f))
  
  src <- basename(map_f) %>% str_remove('_.+')
  spp_map <- read_csv(map_f, show_col_types = FALSE) 

  ### filter for valid cells - AM prob ≥ 0.50, IUCN presence != 5
  prob_thresh <- 0.50
  if(src == 'am') {
    spp_map <- spp_map %>% filter(prob >= prob_thresh)
  } else {
    spp_map <- spp_map %>% filter(presence != 5)
  }
  
  return(spp_map %>% select(cell_id))
}

out_fstem <- here_anx('stressors/max_temp/%s_spp_max_temp_%s.csv')
# zxcv <- list.files(dirname(out_fstem), pattern = 'iucn_', full.names = TRUE)
# list.files(dirname(out_fstem)) %>% n_distinct()

ptm <- proc.time()

for(i in 1:nrow(gps_df)) {
  ### i <- 1
  ### i <- 753
  spp_envs <- T_env_to_process %>%
    filter(gp == i)
  
  tmax_a <- spp_envs$T_max_a[1] ### should all be identical
  tmax_p <- spp_envs$T_max_p[1] ### should all be identical
  
  out_fs <- sprintf(out_fstem, spp_envs$src, spp_envs$id)
  if(all(file.exists(out_fs))) {
    # message(i, ' of ', nrow(gps_df), ': file exists for all ', nrow(spp_envs),
    #         ' species with thermal niche from ', tmax_p, ' to ', tmax_a, '°C...')
    next()
  } else {
    spp_envs <- spp_envs[!file.exists(out_fs), ]
  }
  
  ### not all files exist for this group:
  message(i, ' of ', nrow(gps_df), ': Processing ', nrow(spp_envs), 
          ' species with thermal niche from ', tmax_p, ' to ', tmax_a, '°C...')
  
  t_prs_i_rast <- calc_T_max_prs(abs = tmax_a, pref = tmax_p)
  ### plot(t_prs_i_rast, axes = FALSE, main = sprintf('Therm prs for range %s to %s', tmax_p, tmax_a))
  
  t_prs_i_df <- data.frame(therm_prs = values(t_prs_i_rast),
                           cell_id   = 1:ncell(t_prs_i_rast)) %>%
    filter(!is.na(therm_prs)) %>%
    mutate(therm_prs = round(therm_prs, 5)) ### drop unreasonable precision
  
  tmp <- parallel::mclapply(1:nrow(spp_envs), mc.cores = 16, 
           FUN = function(j) {
    # for(j in 1:nrow(spp_envs)) { 
      ### j <- 2
      ### j <- 119
      ### get species range map
  
      spp_env <- spp_envs[j, ]
      spp <- spp_env$sciname
      src <- spp_env$src
      id  <- spp_env$id
  
      out_f <- sprintf(out_fstem, src, id)
      
      if(file.exists(out_f)) {
        # message('File ', basename(out_f), ' exists... skipping!')
        return(NULL)
      }
      
      ### file doesn't exist; process it!
      message('Processing ', basename(out_f), ' (', j, ' of ', nrow(spp_envs),
        ' in group ', i, ' of ', nrow(gps_df), ')...')
          
      map_stem <- here_anx('spp_maps_mol/%s_spp_mol_%s.csv')
      map_f <- sprintf(map_stem, src, id)
      spp_map <- get_spp_map(map_f)
      
      ### inner join with temp pressure dataframe
      spp_temp_df <- spp_map %>%
        oharac::dt_join(t_prs_i_df, by = 'cell_id', type = 'inner') %>%
        distinct()
      
      ### to check results:
      # spp_range <- map_to_mol(spp_map, by = 'cell_id', which = 'cell_id')
      # plot(spp_range, col = 'black')
      # spp_str <- map_to_mol(spp_temp_df, by = 'cell_id', which = 'therm_prs')
      # plot(spp_str)
      
      ### write out to file
      write_csv(spp_temp_df, out_f)
    })
}

ptm <- proc.time() - ptm

Elapsed time for processing: 3.551 seconds.

3.3.2 Plot example

ex_spp <- c('acropora abrotanoides', 'dermochelys coriacea', 'eustomias achirus', 
            'thunnus albacares', 'prionace glauca')

land_poly <- read_sf(here('_spatial/ne_10m_land/ne_10m_land_no_casp.shp')) %>%
  st_transform(crs = '+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs')
spp_envs <- T_env_rounded %>%
  filter(sciname %in% ex_spp)
for(i in 1:nrow(spp_envs)) {
  ### i <- 1
  spp <- spp_envs$sciname[i]
  src <- spp_envs$src[i]
  id  <- spp_envs$id[i]
  id_lbl <- ifelse(src == 'iucn', paste('iucn', id), 'aquamaps')
  tmax_a <- spp_envs$T_max_a[i]
  tmax_p <- spp_envs$T_max_p[i]
  
  cells <- read_csv(sprintf(out_fstem, src, id)) %>%
    distinct()
  r <- map_to_mol(cells, which = 'therm_prs')
  map_df <- as.data.frame(r, xy = TRUE) %>%
    filter(!is.na(therm_prs))
  
  p <- ggplot() +
    geom_raster(data = map_df, aes(x, y, fill = therm_prs)) +
    scale_fill_distiller(palette = 'RdYlGn') +
    geom_sf(data = land_poly, size = .1, color = 'grey50', fill = 'grey90') +
    theme_void() +
    labs(fill = 'Pressure',
         title = sprintf('%s (%s): T_max_p = %s°C, T_max_a = %s°C', spp, id_lbl, tmax_p, tmax_a))
    
  print(p)
}